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Numerical  treatment  of  the  elliptic  boundary  value  problem  with  nonsmboih" 
solution  by  the  finite  element  method  is  discussed.  The  nonsmoothness  could 
have  its  origin  in  the  unsmooth  boundary  or  the  differential  equation.  This 
paper,  which  is  a  survey  of  the  recent  results,  elaborates  among  others  on 
the  method  of  auxiliary  mapping,  the  partition  of  unity  finite  element  method 
and  the  hp  version  of  FEM  in  three- dimensions.  Numerical  examples  illustrate 
mathematical  results. 


1  Introduction 

Finite  element  methods  face  significant  problems  if  the  exact  solution  of  the  solved 
problem  is  not  sufficiently  smooth. 

The  unsmoothness  could  have  very  different  character. 

a)  Let  us  consider,  eis  the  model  problem,  the  boundary  value  problem  for  the 
Laplace  or  elasticity  equation  on  the  domain  C  R^.  If  the  boundary  of  the 
domain  has  corner,  located  with  the  internal  angle  uj  in  the  origin,  then  the 
solution  u,  in  the  neighborhood  of  the  origin,  has  essentially  the  form  u  = 
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r^4>{6)  where  {r,9)  are  the  polar  coordinates,  ^  >  0  and  (f>{0)  is  an  analytic 
function  in  9.  For  the  Laplace  equation,  the  mixed  boundary  conditions  and 
internal  angle  u;  =  2?r,  we  have  ^  =  1/4,  For  the  elasticity  equation,  mixed 
boundary  conditions  and  a  particular  internal  angle,  the  exponent  /3  could  be 
very  small,  e.g.  see  [1].  The  same  occurs  when  operators  with  piecewise  smooth 
coefficients  are  considered,  e.g.  see  [2].  The  solution  which  essentially  has  the 
form  u  =  r^<f>{9)  belongs  to  the  Besov  space  for  1  <  A:  <  yS,  but  not 

for  A:  >  1  +  ^.  Also,  u  €  for  A:  <  1  +  /3,  but  u  ^  for  A:  >  1  + 

Furthermore,  the  solution  belongs  to  countably  normed  spaces  introduced  in 
[3],  [4],  and  [5].  Particularly,  we  have  for  any  a  =  (ai,  a2) 

|D"u|  (1.1) 

with  7  >  0,  |a|  >  1  and  C  independent  of  a. 

We  have  shown  in  [3],  [4],  and  [5]  that  the  solution  of  elliptic  problems  with 
piecewise  analytic  data  belongs  to  these  spaces.  In  [6],  [7],  and  [8],  we  have 
proven  that  if  the  exact  solution  belongs  to  this  space  then  the  hp  version  of 
the  finite  element  method  converges  exponentially  with 

||e||s  <  (1.2) 

where  N  is  the  number  of  degrees  of  freedom  and  p  >  0  depends  on  7  in  (1.1); 
particularly,  p  ^  0  as  7  — >  0.  In  [9]  we  have  shown  that  in  one-dimension, 
p  — ♦  0  as  7  — +  0  for  any  hp  version,  more  precisely,  ||e||£;  >  (for  more 

see  §  2.1). 

If  or  7  is  small,  then  we  will  speak  about  a  strong  singularity,  which  should 
be  distinguished  from  the  case  when  or  7  is  not  extremely  small.  We  mention 
here  only  the  two-dimensional  problem,  although  the  situation  is  similar  in 
three-dimensions,  too.  (See  §4  of  this  paper). 

b)  The  second  type  of  unsmooth  solutions  has  completely  different  character.  Con¬ 
sider  the  model  problem  solving  the  equation 


—  Au  -t-  Cu  —  f 


(1.3) 


with  C  being  large  (in  absolute  value)  either  positive  or  negative  (for  example, 
if  considering  Helmholtz  problem)  or 


E 


d  ,  \du 


=  / 


(1.4) 


where  aij(x)  are  rough  functions.  Then  also,  if  /  is  smooth,  the  solution  of  (1.3) 
has  a  boundary  layer  character  (C  >•  1)  or  is  highly  oscillatory  when  C  •<  —  1. 
The  solution  of  (1.3)  or  (1.4)  is  not  smooth,  but  has  very  different  character 
when  compared  with  the  one  described  in  a. 


In  this  paper,  we  will  address  methods  for  solving  problems  when  the  solution  is 
singular  in  the  sense  mentioned  above. 
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2  The  method  of  auxiliary  mapping 

Before  we  address  the  method  of  the  auxiliary  mapping  (MAM),  we  will  discuss  in 
detail  the  one-dimensional  finite  element  method  based  on  polynomial  approximation. 
This  will  lead  to  an  insight  into  difficulties  when  the  standard  finite  element  method 
is  used. 


2.1  The  hp  version  of  the  finite  element  method  in  one- 
dimension. 

Let  us  consider  the  problem 


-u"p  =  f,  x€/=(0,l), 

«/3(0)  =  0,  (2.1) 

up{l)  =  1 

with  the  exact  solution 

Up  =  x^,  1/2.  (2.2) 

This  case  is  analogous  to  the  two-dimensional  case  when  u  =  r°d)(0)  with  a  >  0  and 

/3  =  l/2  +  a. 

Consider  the  finite  element  method  on  the  mesh  A„, 

H  ^  ^  ^  ^  Xji  =1, 

with  Ij  =  (iCj-1,  Xj),  hj  —  Xj  —  Xj_i,  y  =  1,  . . . ,  n  and  polynomial  shape  functions  of 
degrees  pj  >  1  on  Ij,  j  =  1,  . . . ,  n,  with  P„  =  (p,, . . .  ,p„).  Let 

5(A„,P„)  =  {ue  H\l),u{0)  =  0,u(l)  =  1}. 

Denoting  A^(5'(A„,  Pn))  =  dim  5(A„,  P„)  =  —1  -f-  p,-,  N  is  obviously  the  number 

of  degrees  of  freedom. 

Let  be  the  classical  finite  element  solution  of  (2.1),  (2.2).  Furthermore, 

let 

®s(A„,p„)  —  up-  (2.3) 

and 

“^^(A"’-^")  =  ll®SA„.p„)lk=  .  (2.4) 

We  will  now  address  the  question  as  to  what  can  be  said  about  Pp(An,  P„)  as 
a  function  of  (A„,P„).  The  first  essential  question  is  about  the  lower  bound  of 
Ep{An,P„)  with  dim5(A„,P„)  =  N.  In  [9],  we  have  proven  among  others 
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Theorem  2.1  Let 


e^{N)  = 


inf 

dimS(A„,Pn)=N 


■£'/j(An,  Pn)i 


then 


where 


ep{N)  >  C{/3) 


y/(0-l/2)N 
9f _ 

y/W^  ' 


and  C{0)  is  a  constant  independent  of  N. 


□ 


Theorem  2.1  shows  that  when  ^  — 1/2  is  very  small,  i.e.,  if  the  solution  is  strongly 
singular,  then  no  finite  element  method  using  polynomial  shape  functions  can  be  effi¬ 
cient.  Let  us  show  now  that  we  can  construct  a  sequence  (A„,  P„)  so  that  iJ/3(An,  jPn) 
is  essentially  the  same  as  ep{N).  To  this  end,  for  0  <  g  <  1,  define 


A„(5)  :  0  =  Io  <  <  •  •  •  <  a:„ 


with 

Xi  =  q^~\i  =  1, . . . , n 
and  for  0  <  s  let  Pn{s)  =  (pi, . . .  ,pn)  with 

p,-  =  [l  +  s(z-l)],i  =  l,...,n 

where  by  [a],  we  denote  the  closest  integer  to  a.  Denote  by  S{s,q,n)  the  associated 
finite  element  space  for  An{q)  and  Pn{s)  and  let  N{s,  q,  n)  =  dimS{s,  q,  n).  As  before, 
we  denote  by  finite  element  solution  and  by  its  error.  Then  we 

have  proven  in  [9J: 


Theorem  2.2  If 

1.  s  >  So,  then 

(2.5) 

2.  s  <  So,  then 

(2.6) 

S.  s  =  So,  then 

(2.7) 

with 

r  _  1  “  V? 

1  +  v^ 

(2.8) 

and 

.„  =  («- 1/2)!^. 

(2.9) 
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Furthermore,  if 

gopt  =  (V5-1)'  (2.10) 

5opt  =  2a  —  1,  (2.11) 

then  _ 

^  (2.12) 

The  constants  C  in  (2.5),  (2.6),  (2.7),  and  (2.12)  are  independent  of  N.  □ 

Remark  2.1  We  have  proven  in  [9]  that  instead  of  <  used  in  Theorem  2.2,  the 
equivalency  holds,  i.e.,  that  the  lower  bound  heis  the  same  form  as  the  upper  bound. 

In  adaptive  procedures,  often  g  =  1/2  is  used.  Then  for  optimal  distribution  of 
the  degrees  given  by  (2.9), 

(2.13) 

Remark  2.2  We  have  gopt  =  0.1715.  It  is  always  better  slightly  to  overrefine  the 
mesh  and  the  value  g  =  0.15  is  recommended. 


Comparing  (2.12)  with  (2.4),  we  see  that  the  mesh  A„(g)  and  degrees  distribu¬ 
tion  Pn{s)  leads  essentially  to  the  best  possible  error.  To  obtain  the  error  which  is 
comparable  with  the  best  error  (2.4),  it  is  essential  that  degrees  of  elements  are  not 
uniform,  they  are  low  when  the  elements  are  small  and  high  when  the  elements  are 
large.  The  question  arises  what  is  the  error  for  uniform  p.  In  [9],  we  have  proven 


Theorem  2.3  Let 
1.  If  s  >  So,  then 


2.  If  s  <  So,  then 


3.  If  s  =  So,  then 


where 


Pn{s)  =  {pi,.--,Pn),Pi  =  M 


rVlN 


-\/(P-lli)Ny/h:iq\iiT 


<  C(M - ^ 

Ing 


Inr 


So  =  (^-1/2) 

l-v/9 
r  = - - — 

l  +  ^/9 

a  =  min(2/3  —  1,(3). 


(2.14) 

(2.15) 

(2.16) 

(2.17) 

(2.18) 
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Furthermore,  if 


Q  =  Qopt  = 

s  =  5opt  ^  2/3  1, 


(2.19) 

(2.20) 


/{l3-l/2)N/2 


Vn’ 


_  g(<9)  -1/2)N 

~  Vn^ 


(2.21) 


Comparing  Theorem  2.2  and  Theorem  2.3,  we  clearly  see  that  uniform  degrees 
distribution  decreases  the  rate  of  convergence. 

In  [9],  we  analyzed  in  great  detail  the  one-dimensional  case.  From  this  analysis,  it 
follows,  that  for  the  optimal  hp  version,  the  elemental  errors  should  be  asymptotically 
equal,  and  the  degree  in  the  first  element  has  to  be  the  smallest  one.  If  it  is  not,  then 
the  mesh  is  underrefined.  If  the  singularity  is  strong,  the  optimal  mesh  is  so  strongly 
refined,  that  implementational  difficulties  (round-offs)  could  occur. 

As  we  stated  in  Theorem  2.2  and  Theorem  2.3  w  instead  of  <  could  be  used. 

Let  us  illustrate  the  effectiveness  of  the  estimates  in  Theorem  2.3.  Consider  the 
case  P  =  0.75,  q  =  0.15  and  s  =  2^  —  1  =  0.5  and  uniform  p.  This  case  is  analogous 
to  the  case  of  Laplace  equation  mentioned  in  the  introduction  when  7  =  1/4.  Based 
on  Theorem  2.3,  we  have 


Hr.  <  r  ^  --0.5778jV^/^ 

IPS(4,0.1S,n)ll£?  ^  ^^0.375® 

In  Table  2.1,  we  report  for  n,  p,  N  the  relative  error 


(2.22) 


=  l|e^.,„.,5„)ll£/|l«Ml|  ■  100 


and  the  constant 

D{N)  =  (2.23) 

which  illustrates  the  effectiveness  of  (2.22).  Based  on  the  theory  presented  in  [9], 
D{N)  should  be  uniformly  bounded  from  below  and  above.  Because  s  =  1/2,  we 
have  for  odd  n,  pn  =  integer -f  1/2  and  hence,  we  report  the  errors  for  both  p  integers 
closest  to  Sn-  Finally,  we  report  the  ratio  of  /imin/^max-  We  see  that  in  fact,  the 
formula  (2.22)  is  very  accurate 
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n 

P 

N 

e% 

D 

^min/  ^max 

1 

2 

1 

23.81 

0.424 

1. 

2 

1 

1 

23.10 

0.412 

0.18 

2 

2 

3 

15.10 

0.620 

0.18 

3 

1 

2 

17.60 

0.516 

0.26  (-1) 

3 

2 

5 

9.842 

0.655 

0.26  (-1) 

4 

2 

7 

6.786 

0.650 

0.40  (-2) 

5 

2 

9 

5.140 

0.663 

0.59  (-3) 

5 

3 

14 

3.175 

0.742 

0.59  (-3) 

6 

3 

17 

2.206 

0.691 

0.89  (-4) 

7 

3 

20 

1.688 

0.688 

0.13  (-4) 

7 

4 

27 

1.074 

0.744 

0.13  (-4) 

8 

4 

31 

0.6865 

0.621 

0.20  (-5) 

9 

4 

35 

0.5572 

0.645 

0.30  (-6) 

9 

5 

44 

0.3223 

0.616 

0.30  (-6) 

10 

5 

49 

0.2421 

0.595 

0.45  (-7) 

Table  2.1.  The  performance  of  the  hp  version  of  the  finite  element  method. 

We  have  shown  that  Theorem  2.1,  Theorem  2.2,  and  Theorem  2.3  completely  and  very 
precisely  characterize  the  performance  of  the  hp  version.  We  see  that  for  ^0  —  1/2  very 
small,  the  hp  version  using  polynomial  shape  function  cannot  be  effective.  Although, 
the  lower  bound  is  available  only  in  the  one-dimensional  case,  it  is  obvious  that  we  can 
expect  in  two-  (and  three-)  dimensions  that  when  7  is  very  small,  any  finite  element 
method  based  on  polynomial  shape  functions  will  converge  very  slowly. 


2.2  The  2-dimensional  model  problem 

Consider  the  elasticity  problem  (isotropic  material,  E  =  1000,  v  —  0.3)  without  body 
forces  on  the  domain 

n  =  {Xl,  X2  I  l®l|  <  2,  |X2|  <  2}  \  {Xl,  X2  I  0  <  Xi  <  2,  X2  =  0} 

shown  in  Figure  2.1  and  with  the  following  boundary  conditions 

Un  =  Ut  —  O  (fixed)  on  Fj  U  r2  (2.24a) 

Tn  =  10,  Tt  =2  (prescribed  tractions)  on  Fs  (2.24b) 

r„  =  Tt  =  0  (free)  on  dO,  \  (Fi  U  F2  U  F5)  (2.24c) 

The  solution  hcis  major  singularities  in  the  neighborhood  of  the  origin  Pi  =  (0, 0), 
where  7  «  0.3  and  at  the  point  P2  =  (2,2),  where  7  «  0.7.  Hence,  we  have  the 
roughly  analogous  case  as  discussed  in  §  2.1.  In  Figure  2.2(a)  and  Figure  2.2(b),  we 
show  the  mesh  7  with  22  elements  and  mesh  II  with  48  elements.  The  geometric 
factor  is  0.15  (figure  is  not  in  scale),  so  that  the  ratio  between  the  size  of  the  smallest 
and  largest  element  in  the  mesh  II  is  0.843  x  10“^. 
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Figure  2.1:  The  scheme  of  the  model  problem. 


>|  I"  0.15x0.5 


(^)  (b) 

Figure  2.2:  The  mesh  I  (22  elements)  and  mesh  II  (48  elements). 


In  Table  2.2,  we  give  relative  error  in  the  energy  norm  as  a  function  of  the  degree 
p  (uniform)  of  the  elements.  The  error  was  computed  from  the  strain  energy  of  the 
finite  element  solution  and  the  exact  solution  computed  by  extrapolation. 


Mesh  I 

Mes 

\ill 

k 

DOF 

e% 

DOF 

e% 

1 

38 

64.60 

92 

46.52 

2 

120 

49.04 

280 

21.21 

3 

226 

42.26 

488 

15.46 

4 

376 

37.59 

792 

12.12 

5 

570 

34.33 

1192 

10.70 

6 

808 

31.87 

1688 

9.82 

7 

1090 

29.88 

2280 

9.14 

8 

1416 

28.23 

2960 

8.61 

Table  2.2.  The  relative  error  e  in  the  energy  norm  in  %. 

From  Table  2.2,  we  see  an  algebraic  rate  and  that  for  p  >  5,  the  mesh  II  is  under¬ 
refined. 


2.3  The  method  of  auxiliary  mapping 

In  the  neighborhood  of  the  singular  point,  for  example,  the  origin,  the  elements  have 
two  straight  lines  and  one  circular  arc.  The  shape  functions  are  mapped  polynomials 
when  the  blending  mapping  is  used  so  that  the  mapping  of  the  side  {AB)  of  the 
standard  element  on  the  circular  axe  is  lineax  in  the  arc  length  (for  more,  see  [10]). 
The  scheme  is  shown  in  Figure  2.2 


Figure  2.3:  The  elements. 

Let  T  be  the  element  with  one  circular  side  as  shown  in  Figure  2.3. 

Consider  the  conformal  map  z  =  (p  which  mapps  T  onto  T.  Then  we  have  the 
following  lemma: 
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Figure  2.4:  Mapped  element. 

Lemma  2.4  Let  u  be  defined  on  T  and  U  on  T  and  let  u{z)  =  U{C).  Then  for 
i  =  1,2 


xoith  Cl  and  C2  independent  ofu.  □ 

Let  {r,6)  be  the  polar  coordinates  in  the  plane  xi,  X2  and  (p,  V*)  in  the  plane  ^1, 
^2-  Then,  in  the  special  case  u  =  we  get  U  =  p'*°‘ with  for  7  >  0.  Hence 

U  has  smaller  singularity  than  u.  Using  the  pull-back  polynomials  from  the  standard 
element,  we  see  that  U  on.  T  can  be  approximated  in  if^(T)-seminorm  much  better 
them  u  on  T.  In  general,  U  is  smoother  than  u  when  u  is  singular  due  to  the  comer 
of  the  boundary.  Using  Lemma  2.4,  we  can  approximate  u  well  on  T  by  pull-back 
polynomials  if,  in  addition,  the  conformal  map  is  used.  To  obtain  conforming  elements 
with  common  circular  side,  the  blending  mapping  preserving  the  length  as  explained 
above,  is  essential. 

Implementation  of  this  method,  which  we  call  Method  of  Auxiliary  Mapping 
(MAM),  is  very  easy.  This  method  is  very  effective  when  the  solution  is  very  sin¬ 
gular  due  to  corner  or  interface.  In  Table  2.3,  we  show  the  error  in  the  energy  norm 
for  the  MAM  on  the  mesh  I  using  7  =  6  for  the  elements  in  the  origin  Pi  =  (0, 0) 
and  7  =  2  in  for  the  elements  in  the  point  =  (2,2). 
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p 

DOF 

e% 

1 

38 

44.13 

2 

120 

15.81 

3 

226 

9.69 

4 

376 

4.22 

5 

570 

1.80 

6 

808 

1.10 

7 

1090 

0.65 

8 

1416 

0.37 

Table  2.3.  The  accuracy  of  MAM. 


The  MAM  method  on  mesh  /  has  obviously  the  same  DOF  as  the  standard  p-method 
ajid  the  cost  is  also  exactly  the  same.  The  effectiveness  of  the  MAM  method  is  obvious 
by  comparing  Table  2.2  and  Table  2.3.  It  is  necessaxy  to  underline  that  although  the 
method  is  based  on  conformal  mapping,  it  has  nothing  to  do  with  solving  Laplace 
equation.  The  mapping  is  used  for  a  smoothening  and  it  is  essential  that  the  H^- 
seminorm  is  preserved  when  going  from  the  mapped  element  to  the  original  one. 

The  MAM  was  discussed  and  analyzed  in  [11],  [12],  [13]  in  the  two-dimensions 
and  in  [14]  in  three-dimensions.  It  was  shown  that  the  method  is  not  too  sensitive 
to  the  selection  of  the  smoothening  parameter.  The  MAM  is  very  ecisy  to  implement 
as  p  version  or  hp  version  in  the  frame  of  standard  codes.  It  avoids  the  problem  of 
slow  convergence  of  the  classical  hp  version  when  the  solution  is  very  singular  due  to 
corners  of  the  domain  or  interfaces. 


3  The  PUFEM— The  Partition  of  Unity  Finite 
Element  Method 

3.1  Introduction 

Let  us  first  address  the  main  idea  of  the  proof  of  the  convergence  of  the  classical  p 
version  of  FEM.  We  construct  a  function  in  the  finite  element  space  which  approxi¬ 
mates  well  the  exact  solution  in  the  energy  norm.  Then  the  finite  element  solution  has 
the  error  which  is  majorized  by  the  approximations  error  of  the  constructed  function. 
The  construction  proceeds  as  follows  (see  [7],  [8],  [15]): 

(a)  Given  the  partition  of  the  domain  ft  into  elements  T,  an  approximation  by 
polynomials  of  degree  p  on  every  single  element  is  constructed.  Because  of 
the  completeness  of  the  polynomials,  the  error  of  approximation  can  be  made 
arbitrarily  small  by  selecting  sufficiently  high  degree  p. 

(b)  Using  the  fact  that  the  exact  solution  belongs  to  JT^(ft)  (the  energy  space), 
the  discontinuity  of  the  approximation  along  the  element  boundary  can  be  esti¬ 
mated  (in  H^I^{dT))  and  a  correction  is  used  so  that  the  constructed  function 
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is  continuous.  A  theorem  ([15],  [16])  on  the  polynomial  extension  from  dT  to 
T  is  utilized. 

Let  us  assume  now  that  we  solve,  for  example,  the  Laplace  problem,  i.e.,  we  know 
a-priori  that  the  solution  is  harmonic.  Then  realizing  that  the  harmonic  polynomials 
axe  complete  in  the  space  of  harmonic  functions,  we  can  approximate  well  the  exact 
solution  on  every  element  by  harmonic  polynomials  only.  Hence,  there  is  no  difference 
when  compared  with  the  step  (a)  above. 

Nevertheless,  the  part  (b)  now  creates  essential  difficulty  so  that  the  idea  of  ap¬ 
proximation  by  harmonic  polynomials  cannot  be  used.  This  is  especially  important 
when  a  still  more  general  complete  set  of  approximation  functions  is  used.  For  ex¬ 
ample,  Bergmann  [17]  an  Vekua  [18]  construct  an  analog  of  harmonic  polynomials 
satisfying  homogeneous  equations  of  second  order  with  analytic  coefficients.  This  can 
be  used  for  example,  when  equation  (1.3)  with  /  =  0  is  considered. 

The  major  step  is  to  construct  a  continuous  function  from  the  piecewise  continuous 
(on  patches)  functions.  This  will  be  made  by  a  partition  of  unity  approach.  The  idea 
was  used  and  theoretically  analyzed  in  [19],  [20],  [21],  and  [22].  This  is  the  PUFEM 
method  presented  here.  For  detailed  description  of  the  PUFEM  method,  see  [20], 
[21],  and  [22]. 


3.2  The  PUFEM  method 

Let  us  describe  the  major  ingredients  of  the  PUFEM  method.  The  critical  notion  is 
here  the  notion  of  (M,  Coo,  Cq)  partitions  of  unity. 

Definition  3,1  Let  fi  C  R"  an  open  set,  {fl,}  be  an  open  cover  of  Cl  satisfying  a 
pointwise  overlap  condition: 


3M  6  N,  Vx  6  n,  card  {i  |  x  €  }  <  M. 


Let  be  a  Lipschitz  partition  of  unity  subordinate  to  the  cover  {fl,}  satisfying 


supp(^i  C  fij,  Vi, 
^  =  1  on  fi. 


||^»||l<»(r")  <  Coo, 


Cg 

diam  Cl,  ’ 


where  Coo  and  Co  are  two  constants.  Then  {<f>i}  is  called  a  (M,  Coo,Cg)  partition  of 
unity  subordinate  to  cover  {0,}.  The  partition  of  unity  is  said  to  be  of  degree 
ru  €  No  C  C'"(R’").  The  covering  sets  {fl,  }  are  called  patches.  □ 
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Definition  3.2  Let  {fit}  be  an  open  cover  0/  fi  C  R”  and  let  {</>{}  be  a  (M,  Coo,  C'g) 
partition  of  unity  subordinate  to  the  cover  {fi,}.  Let  Vi  C  H^{Cli  D  fi)  be  given.  Then 
space 

V  :='£<f>iVi  =  f^^<f>iVi  1  €  Vij  C  H\n)  (3.1) 

is  called  the  PUFEM  space.  The  PUFEM  space  V  is  said  to  be  of  degree  m  €  No  z/ 
V  C  ^"(fi).  The  spaces  Vi  are  referred  to  as  the  local  approximation  spaces.  □ 

Let  us  now  mention  the  basic  theorem  [20],  [21],  [22]: 

Theorem  3.1  Let  fi  C  R”  6e  given.  Let  {fi,},  {<^{}  and  {1^}  be  as  in  Definition  3.1 
and  Definition  3.2.  Let  u  €  be  the  function  to  be  approximated.  Assume  that 

the  local  approximation  spaces  Vi  have  the  following  approximation  properties:  On 
each  patch  fi,-  H  fi,  u  can  be  approximated  by  a  function  u,-  €  Vi  such  that 


Then  the  function 
satisfies 


Ik  -  v,-|U2(ninn)  <ei{i), 

||V(u  —  uOlU^Cn.-nn)  <  £2(0- 

Uaj>  =  '£<l>iVi£V  cH\D) 


(3.2) 

(3.3) 

(3.4) 


||w  -  «ap|U2(n)  <  VmC^ 

l|V(.  -  (e  eHi)  +  4(.V 


1/2 


(3.5) 

(3.6) 

□ 


The  PUFEM  can  be  understood  as  the  h,  p,  or  hp  version.  Consider  for  example 
the  p  version.  Let  {fi,-}  be  the  patches  covering  fi  and  assume  that  the  exact  solution 
u  is  harmonic.  Assume  that  the  spaces  V  are  the  spaces  of  harmonic  polynomials  of 
degree  p.  Then 


£i(*)  <  Cdiam(ni)p  '•|lu||H»(l!,nn) 
£2(1)  <  Cp"*‘''*||u||Hi(ninQ) 


(3.7) 

(3.8) 


and  the  error  estimate  from  Theorem  3.1  takes  the  form 


l|V(u  -  <  M^CiCa  +  (3.9) 

We  see  that  the  PUFEM  method  allows  to  do  what  was  mentioned  in  §  3.1. 

In  [22],  more  details  of  PUFEM  method  including  the  a-posterion  error  estimation 
are  given. 
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3.3  A  numerical  example 

Let  us  consider  the  Helmholz  problem  on  unit  square 


—  Au  —  k^u  =0  on  n  =  (0, 1)  X  (0, 1), 
du 

— — h  iku  —  g  on  dCl, 


(3.10) 

(3.11) 


where  g  is  chosen  such  that  the  exact  solution  u  is  a  plane  wave  of  the  form 

u(a;)  =  ^  . 
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(3.12) 


The  following  types  of  local  approximations  spaces  were  analyzed  in  [21].  The 
first  type  are  “generalized  harmonic  polynomials”  of  Bergmann-Vekua  type.  Then 


V’ip)  =  span  {e*'"V„(J:,  r)  |  n  =  0, . . .  ,p} 


(3.13) 


where  the  functions  Jn  are  Bessel  functions  of  the  first  kind.  The  second  type  are 
systems  of  plane  waves  given  by 

W{p)  =  span  |  =  —jj  =  0, . . .  ,p  -  1}  (3.14) 

P 

with  p  =  4n  +  2,  n  an  integer,  so  that  u{x)  given  by  (3.12)  does  not  belong  to  W{p) 
for  any  p.  It  can  be  shown  that  functions  from  V'^^k)  and  W{p)  can  approximate 
well  any  function  satisfying  (3.10).  For  more,  see  [2lj.  If  u  is  solution  of  (3.10),  and 
O  C  n,  then 

11“  “  “fllif®  -  C(7,fi)e-"' 

““  “  “''*«■<“>  - 

holds  for  any  7  >  0. 

Assume  now  that  on  fi,  we  have  a  square  mesh  with  size  h  =  Ijn  and  use  the 
partition  of  unity  created  by  the  standard  bilinear  pyramid  shape  functions. 

We  now  compare  the  effectiveness  of  the  PUFEM  method  with  some  other  ones. 
Let  us  discuss,  as  example,  the  Ccise  k  =  100  and  be  interested  in  the  error  measured 
in  the  L2  norm.  We  will  compare  the  PUFEM  based  on  the  space  W{p)  with  the 
usual  Gaberkin  method  (FEM),  the  generalized  least  squares  finite  element  method 
(GLSFEM  [23])  and  the  quasi-stabilized  finite  element  method  (QSFEM  [24]).  The 
FEM,  GLSFEM  and  QSFEM  are  based  on  piecewise  linear  functions  on  uniform 
mesh  and  they  differ  in  their  choice  of  bilinear  form.  In  particular,  the  bilinear  form 
of  QSFEM  is  constructed  such  that  the  “pollution”  (see  [24])  is  minimized  and  it  is 
virtually  the  best  method  available  which  is  based  on  piecewise  linear  functions. 

In  Table  3.1,  we  show  the  necessary  DOF  to  obtain  relative  error  e%  in  L2  norm. 
We  report  in  Table  3.1  also  the  DOF  of  the  best  approximant  by  bilinear  functions. 


14 


e% 

Best  approx,  by 
bilin.  funct. 

QSFEM 

GLSFEM 

FEM 

30 

2.045  (+3) 

3.969  (+3) 

2.016  (+4) 

7.734  (+4) 

10 

5.041  (+3) 

1.000  (+4) 

6.150  (+4) 

2.352  (+5) 

5 

3.464  (+3) 

1.960  (+4) 

1.274  (+5) 

4.692  (+5) 

Table  3.1.  DOF  necesary  to  obtain  the  accuracy  e%  in  L2  norm;  k  =  100. 

In  Table  3.2,  we  show  the  DOF  for  PUFEM,  forn  =  4  and  compare  it  with  the  other 
methods 


P 

PUFEM 

Best  approx,  by 
bilin.  funct. 

QSFEM 

FEM 

26 

10.8 

6.50  (+2) 

3.80  (+3) 

7.95  (+3) 

2.08  (+5) 

30 

0.69 

7.50  (+2) 

5.89  (+4) 

1.23  (+5) 

3.23  (+6) 

34 

0.11 

8.50  (+2) 

3.45  (+5) 

7.23  (+5) 

1.90  (+7) 

Table  3.2.  DOF  necesary  to  achieve  various  accuracies  in  L2  for  PUFEM  with 
n  =  4  and  various  other  methods,  k  =  100. 

The  excellent  performance  of  the  PUFEM  method  shown  in  Table  3.2  is  due  to 
the  fact  that  PUFEM  employs  the  general  character  of  the  solution;  while  the  other 
methods  utilize  only  piecewise  linear  functions.  In  Table  3.3,  we  report  the  number 
of  operations  using  band  elimination,  for  PUFEM,  n  =  4. 


P 

e% 

PUFEM 

QSFEM 

FEM 

26 

10.8 

1.76 

(+7) 

6.3  (+7) 

4.3  (+11) 

30 

0.69 

2.71 

(+7) 

1.5  (+10) 

1.01  (+13) 

34 

0.11 

3.94 

(+7) 

5.2  (+11) 

3.6  (+14) 

Table  3.2.  The  number  of  operations  using  band  elimination,  k  =  100,  error  in  L2, 

k  =  100. 

Let  us  underline  that  the  construction  of  the  stiffness  matrix  for  the  PUFEM  based 
on  W{p)  is  relatively  cheap  and  the  cost  is  negligible  when  compared  with  the  cost 
of  the  solver.  The  stiffness  matrix  for  the  PUFEM  based  on  the  space  V"(p)  is  more 
expensive  than  for  W{p).  Nevertheless,  the  space  V'’{p)  is  in  some  sense  optimal  (see 
[22]). 

In  Table  3.4  and  Table  3.5,  we  show  the  results  for  k  =  Z2  and  the  error  measured 
in  the  if^-seminorm.  Table  3.4  shows  the  accuracies  and  number  of  iterations  and 
operations  for  the  iterative  method  proposed  in  [25].  Table  3.5  shows  the  operation 
count  for  PUFEM  method  with  the  band  elimination,  k  —  32  and  n  =  1. 
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FEM 

,'SFEM 

Vdof 

error 

#of 

NOP 

H'’  error 

#of 

NOP 

e% 

iter  at. 

e% 

iterat. 

33 

65.0 

232 

4.51  (+6) 

30.5 

272 

5.29  (-1-6) 

65 

21.7 

434 

3.37  (-1-7) 

14.3 

492 

3.82  (-1-7) 

129 

8.16 

831 

2.68  (-1-8) 

7.02 

953 

2.96  (4-8) 

257 

3.64 

1665 

2.07  (-f9) 

3.48 

1863 

2.31  (-1-9) 

513 

1.72 

3263 

1.62  (-flO)  1 

1.69 

3762 

1.86  (-1-10) 

Table  3.4.  Operation  count  for  solving  the  linear  system,  error  in  /f^-seminorm 

k  =  32. 


P 

error  e% 

NOP 

18 

46. 

1.3  (-1-5) 

22 

6.7 

2.3  (-1-5) 

26 

0.38 

3.8  (-h5) 

30 

0.00025 

5.9  (-1-5) 

Table  3.5.  Operation  count  for  band  elimination  for  PUFEM,  A;  =  32,  n  =  1,  error 

in  -seminorm 

Let  us  underline  once  more  that  the  cost  of  the  stiffness  matrix  is  negligible  in  compar¬ 
ison  with  the  cost  of  the  solver.  The  PUFEM  method  belongs  to  the  family  of  mesh 
free  methods,  see  [26],  [27],  nevertheless,  we  underline  in  this  paper  the  flexibility  of 
PUFEM  which  allows  us  to  employ  the  properties  of  solved  differential  equations. 


4  Th©  hp  version  of  FEIVI  in  three- dimensions 
4.1  The  hp  version 

In  contrast  to  the  two-dimensional  case,  the  character  of  the  singularities  in  the 
neighborhood  of  the  boundary  is  much  more  complex.  In  three-dimensions,  we  have 
to  distinguish  between  the  behavior  in  the  neighborhood  of  the  edges  far  from  the 
vertices,  close  to  the  vertices,  and  in  the  neighborhood  of  the  vertices,  and  in  the 
neighborhood  of  the  vertices  which  is  (conically)  far  from  the  edges.  In  the  two- 
dimensional  case,  only  one  type  exist  in  the  neighborhood  of  the  corner. 

In  [28],  [29],  [30],  and  [31],  the  regularity  of  the  solutions  in  terms  of  countable 
spaces  is  analyzed.  Based  on  these  results,  the  hp  versions  converge  exponentially 

||e||<Ce-^^^'*,7>0.  (4.1) 

In  contrast  to  the  one-dimensional  case,  there  is  no  proof  that  the  exponent  1  /5 
IS  the  optimal  one,  but  we  conjecture  that  it  is.  The  meshes  leading  to  this  rate 
have  elements  with  large  aspect  ratios,  which  increases  as  TV  oo.  These  “needle” 
elements  are  in  the  neighborhoods  of  the  edges  and  reflect  the  fact  that  the  solution 
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is  smooth  along  the  edges  and  unsmooth  in  the  direction  perpendicular  to  the  edge. 
Also,  the  optimal  upper-bound  (and  the  lower-bound)  depending  on  the  strength  of 
the  singularities,  is  not  known,  but  analogously  as  in  the  one- dimensional  case,  it  is 
possible  to  expect  that  7  in  (4.1)  could  be  very  small  in  the  case  when  the  singularity 
is  strong.  So  far,  we  discussed  the  hp  version.  The  p  version  was  analyzed  in  [32]  and 
[33].  See  also,  [34]  for  some  computational  experience. 

4.2  Computational  Examples 

In  this  section,  we  show  briefly  an  example  analyzed  by  the  code  STRIPE  developed 
at  Aeronautical  Research  Institute  of  Sweden.  Consider  the  problem  on  the  domain 
shown  in  Figure  4.1  where  the  boundary  conditions  are 

•  on  the  faces  A-C-E-N,  A-B-E-F,  A-B-C-D,  and  I-J-L-M,  we  have 

(r„r5„r,)  =  o; 

•  on  the  face  G-H-I-J,  we  have  {Tx,Ty,Tz)  =  (—10,0,0); 

•  on  the  face  D-K-LrM-N-C,  we  have  u  =  0,  Ty  =  =  0] 

•  on  the  face  N-M-J-H-F-E,  we  ahve  v  =  1,  Ti  =  =  0; 

•  on  the  face  K-G-I-L,  we  ahve  u  =  0,  Tr  =  T*  =  0;  and 

•  on  the  face  K-D-B-F-H-G,  we  ahve  w  =  0,  Tx  =  Ty  =  0. 


Figure  4.1;  The  domain 


The  singularity  of  the  solution  occurs  in  the  vertex  A,  along  the  edges  A-B,  A-C, 
A-E,  and  the  edge  I-J.  The  singularity  along  the  edge  I-J  is  weaker  than  along 
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A-B^  A-C,  and  A-E.  The  basic  mesh  is  shown  also  in  Figure  4.1.  Around  the  edges, 
the  mesh  is  refined  in  the  analogous  way  as  in  two-dimensions,  with  Mi,  M2  layers 
in  the  neighborhood  of  the  edges  A-B,  A-C,  A-E,  and  I-J.  Figure  4.2  shows  the 
case  with  Mi  =  1  and  M2  =  3.  In  the  neighborhood  of  the  vertex  A,  the  mesh  is  a 


Figure  4.2:  Detail  of  the  Mesh 


complex  one  accommodating  the  meshes  along  the  edges  A-B,  A-C,  and  A-E.  The 
hp  version  together  with  increasing  N  does  increase  the  number  of  layers  around  the 
edges  where  the  singularity  is  located.  The  p  version  fixes  the  mesh  and  increases  the 
degree  uniformly  or  selectively.  In  the  Figure  4.3,  we  show  in  the  scale  log  ||e||  x 
the  error  for  various  p  (uniform)  and  Mi  with  M2  =  1.  As  predicted,  we  see  the  rate 
if  we  combine  properly  Mi  and  p.  It  will  appear  in  the  graph  as  the  straight 
line.  As  seen  from  Figure  4.3,  we  have  to  combine  in  a  proper  way  the  degrees  of 
elements  and  the  number  of  layers.  In  the  range  of  accuracy  under  consideration, 
M2  =  1  is  sufficient.  In  Table  4.1,  we  show  the  proper  combination  which  leads  to 
the  exponential  convergence  with  respect  to  N  =  DOF 


Ml 

P 

0 

6 

2 

5,6,7 

4 

6,7,8 

Table  4.1.  The  combination  of  Mi  x  p  for  the  hp  version. 

Remark  4.1  The  relative  error  in  the  energy  norm  shown  in  Figure  4.3  was  com¬ 
puted  from  the  strain  energies 

II® Hi:  / ®exact 

||w|1e  V  eex  y 


18 


p-VERSION 


Figure  4.3:  The  error  in  the  energy  norm 


where  the  energy  £exact  was  computed  from  Mi  =  4,  M2  =  2  and  p  =  10  and 
extrapolation  so  that  the  error  in  the  given  range  shown  in  Figure  4.3  is  guaranteed. 

□ 

Figure  4.3  shows  also  the  error  for  the  h  version  with  uniform  mesh  and  elements 
of  degree  2.  In  Figure  4.4,  we  show  the  CPU  time  on  one  processor  computation  of  a 
Silicon  Graphics  Challenger  for  the  hp  version  and  the  h  version.  In  the  CPU  time, 
the  total  time  is  included,  i.e.,  the  construction  of  the  stiffness  matrix,  direct  solver 
and  the  postprocessing.  In  Figure  4.4,  we  show  also  CPU  time  for  the  iterative  solver 
written  by  J.  Mandel  (Solver  International,  Inc.)  which  is  a  special  solver  based  on 
PCG  method. 

So  far,  we  have  addressed  the  case  with  uniform  degrees.  We  can  also  use  an 
adaptive  procedure  which  is  more  effective.  In  the  adaptive  procedure,  different 
degrees  of  the  shape  functions  (but  not  the  full  space  of  shape  functions  of  degree  p) 
are  used.  In  Figure  4.5,  we  show  the  relation  between  the  CPU  time  and  achieved 
accuracy  measured  in  the  energy  norm  for  the  adaptive  approach  and  uniform  p 
approach  for  Mi  =  4  and  M2  =  1. 

We  see  clearly  the  effectiveness  of  the  hp  version  in  the  three-dimensional  examples 
for  solutions  with  singularities. 
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